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Abstract. We consider the method of particular solutions for numerically 
^ ' computing eigenvalues and eigenfunctions of the Laplacian on a smooth, bound- 

ed domain Q in M" with either Dirichlet or Neumann boundary conditions. 
fvj , This method constructs approximate eigenvalues E, and approximate eigen- 

functions u that satisfy Au = Eu in f2, but not the exact boundary condition. 
An inclusion bound is then an estimate on the distance of E from the actual 
spectrum of the Laplacian, in terms of (boundary data of) u. We prove op- 
CU , erator norm estimates on certain operators on L^(dfl) constructed from the 

r^) I boundary values of the true eigenfunctions, and show that these estimates lead 

to sharp inclusion bounds in the sense that their scaling with E is optimal. 
This is advantageous for the accurate computation of large eigenvalues. The 
jrt , Dirichlet case can be treated using elementary arguments and will appear in 

[5], while the Neumann case seems to require much more sophisticated tech- 
nology. We include preliminary numerical examples for the Neumann case. 
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1. Introduction 



In this paper we consider Laplace eigenfunctions on a smooth bounded domain 
^^ , Q C K". As is well known, the positive Laplaciary, 

i=l ' 



with domain either H^{U) n Hq{^) or 



H I is self-adjoint. Here and below, d„ denotes the directional derivative with respect 

to the outward unit normal vector at dfl. These are known as the Laplacian 
with Dirichlet, resp. Neumann, boundary conditions and will be denoted Au, 
resp. An- In either case, there is an orthonormal basis of L^{n) consisting of 
real eigenfunctions. We denote an orthonormal basis of Dirichlet eigenfunctions Uj, 
j = 1 ... 00, and an orthonormal basis of Neumann eigenfunctions Vj, j = 1 ... 00, 

A,j , Vj satisfy a 



with Dirichlet/Neumann eigenvalues Ej = 


A|, resp. 


E,= 


M? 


. Thus 


u 


Helmholtz equation 














Auj = 


= EjUj — XjUj, or 


Avj = 


Ejv, 


= 


MjWj, 




with boundary condition 














Uj\dn = 0, or 


dnVjldfl 


= 0. 









Note that our sign convention is opposite to that of [5] 
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Figure 1. Graph of the minimum value of t[u] achievable at each E, 
for the Dirichlet case; see (|l.ip . The smooth domain f2 C R^ is shown in 
Fig. [l] u is restricted to lie in the span of a finite number of numerical 
basis functions satisfying the Helmholtz equation at each E, and the 
minimization uses dense linear algebra [S]. A small value (e.g. at the 
dot shown) implies closeness to a Dirichlet eigenvalue. 



We will denote the spectrum of Ad, resp. A^v by spec^, resp. spec^. Also, we 
will denote the normal derivative dnUj at dil by ipj, and the restriction of Vj to dft 
by Wj. Thus tpj-,Wj are functions on dfl, which we will refer to as the boundary 
traces of eigenfunctions Mj, resp. Vj. 

The Method of Particular Solutions fS", '3' is a numerical method for finding 
eigenvalues and eigenfunctions of the Laplacian on a Euclidean domain. In the case 
of the Dirichlet boundary condition, the method consists of choosing an energy 
(positive real number) E, and looking for the solution u to the Helmholtz equation 
(A — E)u = that comes closest to satisfying the boundary condition, in the sense 
that it minimizes (or approximately minimizes) the L^ norm of the boundary trace 
of M (in L^(9ri)). In practice u is restricted to a sufficiently large numerical subspace 
and the minimization is done via dense linear algebra (a generalized eigenvalue 
or singular value problem [3J [5]). We then think of this minimum L- norm on 
the boundary as a function of E (see Fig. [1]), and numerically try to find the 
(near) zeros of this function as we move along the E-a,xis. Clearly, if we find a 
Helmholtz u with llwHis/gQ), then E is a, Dirichlet eigenvalue, and u is a Dirichlet 
eigenfunction. We would expect, therefore, that if ||M||L2(gQ) is very small, then E 
is close to a Dirichlet eigenvalue, and u is close to (i.e. makes a small angle with) 
the corresponding eigenspace. 

An inclusion bound is a quantitative estimate of this form, taking the form (for 
eigenvalues) 



(1.1) 



d{E, spec jy) <CE"t[u], where t[u] 






for some constant C independent of E and exponent a. This tells us how small 
we must make the 'tension' t[u\ in order to achieve any desired accuracy in our 
numerically computed eigenvalue. We are primarily interested in high energy es- 
timates, i.e. E large, so our goal is to obtain such an estimate which is sharp as 
i? — > oo, that is, with the smallest possible a. (Ideally, for practical applications, 
we also want a constant C that is small and computable, that is, expressed in terms 
of geometric quantities such as the measure, surface measure, inradius, etc., of our 
domain fi. But we will not address that issue here.) 
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This article is a report on completed work on the Dirichlet case |5], and an 
announcement of work in progress on the Neumann case. Full details will appear 
elsewhere. 

2. Dirichlet boundary condition 
We begin by proving that there are upper and lower bounds 

(2.1) C-^\,<U,\\L^(on)<CX, 

where C depends only on fl. These are classical and well-known estimates, but we 
give the proof since it follows from a calculation that we need later on anyway. The 
proof is via a Rellich-type identity, involving the commutator of A with a suitably 
chosen vector field V . The basic calculation is 

(u, [A, V]u) = I (((A - \^)u){Vu) - m(1/(A - X^)u)^ 

(2.2) +/ ({d,,u)(Vu)^u{dAVu))). 

Jan ^ ' 

If u = Uj is a Dirichlet eigenfunction with eigenvalue A^, then three of the terms 
on the RHS vanish, and we obtain 

{u,[A,V]u)^ f id^u){Vu). 
Jon 

If we choose V so that, at the boundary, it is equal to the exterior unit normal, 
then the RHS is precisely HV'jlP- (If not indicated, norms will be assumed to be L^ 
norms.) The left hand side is {u, Qu) where Q is a second order differential operator 
and is O(A^), yielding the upper bound HV'jlP — O(A^). On the other hand, if we 
take V to be the vector field ^^ Xidxi, then [A, V] = 2 A. Then the LHS is exactly 
equal to 2A^, while the RHS is no bigger than (max^o |a::|)|ji/)j|p, yielding the lower 
bound A2 = OiWijjf). This lower bound is due to Rellich [lU] . 

It turns out that there is a very useful generalization of the upper bound in (j2.ip , 
proved recently by the authors, that applies to a whole 0(1) frequency window: 

Theorem 2.1. Let Q C M" be a smooth bounded domain and let tpi be defined as 
above. Then the operator norm of the finite rank operator 

(2.3) Y. Mi''^r) ■.L\dQ)^L\dQ) 

A,G[A,A+1] 

is bounded by CX^ , where C depends only on fl. 

Remarks: 

• This is quite a strong estimate, since by the lower bound of (|2.1|) there is a 
lower bound of the form cX^ on the operator norm of any one term in the 
sum. 

• This is closely related to the phenomenon of 'quasi-orthogonality' of V'i 
and tl^j, when |Ai — A^j is small. Indeed, this estimate implies that when 
|Ai — Ajl < 1, then the inner product (V'i,V'j) is usually small compared 
with A^. See Barnett [2]. 

• This is also closely related to an identity of Backer, Fiirstberger, Schubert 
and Steiner [1]. 
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Theorem 1 2. II is proved as follows: first, we prove the upper bound ||rfnu||L2(an) < 
CA||u||l2(j7-| is valid not just for eigenfunctions, but for approximate eigenf unctions 
u G dom Ad such that 

||(A-A2)«|U.(n) = 0(A). 

In fact, the proof is almost unchanged: we use p.2p again. Now the term (A — 
X'^)u{Vu) is no longer zero, but by assumption (A — A^)u is 0(A) in L^(il), and also 
Vu is 0(A) in L?{Vt)^ so by Cauchy-Schwarz this term is O(A^). We treat the term 
u{V{A — \^)u) similarly, after first integrating the vector field V by parts (which 
produces no boundary term since u vanishes at d^). The rest of the argument runs 
as above. This was also noticed by Xu [E]- Notice that this condition applies in 
particular to a spectral cluster, that is, for u S range -E[a.a+i] (a/^-d)- We then use 
a TT* argument: We define an operator T from range i?[A,A+i] W ^d) to L?{d^) 
by Tu = dnu\Qn- We can express T in terms of the eigenfunctions by 

Tu= ^ {u,Ui)tpi. 

A,G[A,A+1] 

That is, it is just the normal derivative of the element ^ UiUi of this spectral 
subspace. Then, as we have just shown, 

||T|| <CX. 

It follows that TT* : L'^{dn) -> L'^{dn) has operator norm bounded by C^A^ But 
TT* is precisely the operator (|2.3p appearing in the statement of the theorem. 

3. DiRICHLET INCLUSION BOUND 

As mentioned above, a Dirichlet inclusion bound is an estimate on the distance 
from E to spec^ in terms of t[u] defined in (jl.ip of a Helmholtz solution u. A 
classical result along these lines is the Moler-Payne inclusion bound [9 . This says 
that 

d{E, specjy) < CEt[u], (A - E)u = 0, 

where t[u\ is defined in (jl.ip and C depends only on fl. The proof in [9] uses very 
little about the Dirichlet problem in particular. 

Recently, Barnett [3], followed by the authors 5 , improved this bound by a 
factor of \/E: 



Theorem 3.1. There exist constants c,C depending only on fl such that the fol- 
lowing holds. Let u he any nonzero solution of (A — E)u = in C°°{fl), and let 
^*miii be the Helmholtz solution minimizing t[u\. Then 

cVEt[u,nin\ < d{E,specjj) < CVEt[u]. 

Remark 3.2. There is always a Helmholtz solution that minimizes t[u]; see [S]. 

Proof. The result is trivial if E' G specj^A. Suppose that E is not an eigenvalue, 
and consider the map Z{E) that takes / G L^{dfl) to the (unique) solution u of 
the equation 

{A-E)u = 0,u\an^f. 
The u that minimizes t[u] then maximizes ||u||L2m) given |lw||^2(gQ). So 

(3.1) (mintM)-^ = ||Z(£;)|| =^ {mm t[u]y^ = \\A{E)\\, 
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where the operator A{E) : L^{dVl) -^ L^{dn) is defined by A{E) = Z{E)*Z{E). 
We claim that A{E) has the expression [3] 

Here and in the remainder of the article, we take the symbol min to mean the 

u 

minimum over the space of Helmholtz solutions. 

To prove p.2p . we show that Z{E) has the expression 



(3.3) Z{E)! ^ J2 ■ 



E,, 



from which (|3.2I) follows immediately. To express Z{X), suppose / is given and u = 
Z{E)f. We write u = J2 o-iUi as a linear combination of Dirichlet eigenfunctions. 
Then, using the Helmholtz formula and Green's identities, 

a-i = {u,Ui) = / ((Am)mj - u{Aui)) 

= ^ _ ^. / {u{dnUi) - {dnU)ui) = ^_ ^ / fi^i 



which proves 

The lower bound in Theorem 13.11 is easy to prove: we note that A{E) is a sum 
of positive operators in (13.21) . so the operator norm of A{E) is bounded below by 
the operator norm of any one summand. So, using the upper bound in (|2.1|) . 



\ME)\\ > 



i'A^j 



cE 
> 



[E-Ejf - d{E, speco) 



2 ■ 



where Ej is the closest eigenfrequency to E. Since (mint[u])^^ = ||A(£')||, this 
proves the lower bound. 

To prove the upper bound, we use Theorem 12. II We need to show that 

To do this, we break up the sum p.2p into the 'close' eigenfrequencies in the 
interval [A — 1, A + 1] and the rest. The estimate p.4|) for the close eigenfrequencies 
is immediate from Theorem 12. II 

For the far eigenfrequencies, we first treat those that lie in the interval [A/2, 2A]. 
These can be broken up into frequency windows of width 1, which are distance 
1,2,3, etc from the chosen frequency A. Then, in p.2p . the numerator for each 
window have operator norm bounded by E, and the denominator is n^E. Since 
^ n^^ is finite, the contribution from these eigenvalues is 0{1). 

For the the eigenfrequencies not lying in [A/2, 2A], it is not hard to see that the 
contribution is only 0{E~^^^), since we get a factor E-^ in the denominator. So the 
far eigenfrequencies altogether only contribute only 0(1) to the operator norm of 
AiE). 

Finally we observe that 

CE ^^ 

d{E,specjj)'^ 

for every E, since the distance from E to the spectrum can be at most ~ Ve (this 
follows by considering an approximate eigenfunction supported in a ball contained 
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in f2). Thus the contribution from the near eigenvalues dominates, and we see that 
()3.4p is true for all E, completing the proof. D 



Similar reasoning gives a bound for the distance between u and the closest eigen- 
function: 

Theorem 3.3. There is a constant C depending only on Vi, such that the following 
holds. Let E > 1, let Ej be the eigenvalue nearest to E, and let Ek the next 
nearest distinct eigenvalue. Suppose u is a solution of {A + E)u = in C°°{n) with 
ll'"l|L2(t2) = Ij o.iT'd- let Uj be the projection of u onto the Ej eigenspace. Then, 



(3.5) \\u~u,hHn) < C-^*[''] 



\E-Ek 



For the proof, see [5]. 

4. Neumann boundary condition 

We next consider the method of particular solutions for computing Neumann 
eigenvalues and eigenfunctions. The Neumann boundary condition is d„w|ao — 0. 
It seems natural to minimize (cf. ()2.3p ) 

hdiv] = —rr-r. '—^, 

over nontrivial solutions v of {A — E)v = 0, since iid[w] — implies that E is 
a Neumann eigenvalue and v a Neumann eigenfunction. Notice, though, that we 
could equally well minimize the quantity 

r r 1 \\F{dnv)\\L^an) 

for any invertible operator F on L^{dil). It turns out that there is an essentially 
optimal choice of F (which depends on E), which is not the identity. Indeed the 
main point of this article is to determine this optimal F. 

The form of F is suggested by the local Weyl law for boundary values of eigen- 
functions. This law [71 [5] says that the boundary traces of eigenfunctions are, on 
the average, distributed in phase space T*{dfl) according to 

c{l-\v\Y^%M<i} (Dirichlet), 

5(1 - |??n"^/*l{|,,|<i} (Neumann) 

where c, c are constants depending only on dimension. This is in the sense of 
expectation values; that is, for any semiclassical pseudodifferential operator A on 
dO, with principal symbol a{y,r]), {y,ri) e T*dVL, we have 

(4.2) 

^^™ AT i\\ H \~^(V'i, AtiV'j) = c / (1 - \r^\^Y''^l{\,^\<^a{y,T^)dydri, 



A,<A 



— - V {wj,A -iWj) ^c (l-k/P) ^^'^l{\r^\<i}a{y,ri)dydr] 



lim , \"'7,-^ -^-"v/ — ^ ( 
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where Nd,Nn are the Duichlet, resp. Neumann eigenvalue countmg functions. 
Here y S dil, {y,r]) e T*{dn) and JTyp is calculated with respect to the induced 
metric on the boundary. See [8, . 

Here we have adopted the semiclassical scaling, that is the wavevectors at eigen- 
value A^ are scaled by ft, = hj = \J^ or jij^ so that they are rescaled to have length 
1 in M", and therefore length < 1 when restricted to the boundary. 

An intuitive explanation for the difference in the distribution of the ■0^ and 
the Wj is as follows. If we use Fermi normal coordinates (y, r) near 9f2, where 
r is distance to dQ, and if (?7,p) are the dual cotangent coordinates, then the 
symbol of the semiclassical operator h^A — 1 is a{h^A — I) = p^ + |?7p — 1. The 
semiclassical normal derivative ihdn has symbol p, which when restricted to the 
boundary and to the characteristic variety {a{h^A — 1) = 0} is equal to a/1 — jryp. 
Since the Dirichlet expectation value A~ {ipj, A^-iipj) involves the application of 

two semiclassical normal derivatives (one for each factor of ijjj), compared to the 
Neumann expectation value {wj, A -iWj), it is not surprising that the Dirichlet 

distribution in (|4.2[) is 1 — |ryp times the Neumann distribution. We can draw a 
moral from this. 

Moral: semiclassically, in the Neumann case, the boundary trace 
that is analogous to dnU in the Dirichlet case is not w|an, but rather 
(4.3) (1 — /i^Aao)_^' {v\qq), where Ago is the (positive) Laplacian on the 
boundary, m, v are Helmholtz solutions at energy /i^^, and (...)+ de- 
notes the positive part. 

To see what goes wrong with using the naive measure iid of the 'boundary condition 
error', let us attempt to follow the same reasoning as in Section[3] We can certainly 
show that 



(4.4) (miniidM)-' = IE 






3 \" f-J ^ 



L^(an)^L^{dn) 



The problem is that the Wj do not behave as uniformly as the Dirichlet traces ipj] 
we have a lower bound 

(4-5) \\wj\\L^dn)>c, 

but the sharp upper bound is 

(4.6) \\wjh2^gn^<Cpy\ 

(This estimate follows from Tataru [IT,.) The reason why, in Theorem 13.11 we 
were able to get upper and lower bounds on d{E,s'pccjj) of the same order in E 
was that the lower bound on the operator norm of a single term 'ipj{tpj^ ■) was of 
the same order as the upper bound on the sum ^ ■ ipj{il)j^ ■) over a whole spectral 
cluster |A — Ajl < 1. In the Neumann case, using i^ will lead to a gap of at least 
^1/3 = E'l/B between the upper and lower bounds on d{E, spec^). 



Notice, however, that if we take our Moral, (|4.3p . seriously, then we could expect 
to find good upper and lower bounds on the quantity (1 — h?,AQQ)jl_ Wj instead. 
Indeed, this is the case, and we have the following exact analogues of (|2.ip . and 
Theorem [SHI 
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Theorem 4.1. Let f2 C M" he a smooth bounded domain, and let Wj be the re- 
striction to dfl of the jth L^ -normalized Neumann eigenfunction Vj. Then there 
are constants c, C such that 

(i) 11(1 - h'^jA9a)+'^Wj\\L2^Qn) > c, hj = fij^; 
(ii) the operator norm of 
(4.7) Y. {^~h^Aan)T^o {{^-h^^an)T^or), h = fi-\ 

is bounded by C . 

Example 4.2. On the unit disc, Neumann eigenfunctions have the form 

v{r,e) = ce™^J„(/i„,ir), where J^{nn,i) = 0, 
and from (12.21) we derive 



2mL- / U,i-n')\v\' =^ ||(l-Aao/M,')f«^.ll=V2. 

Since zeroes of J,'j are at least n apart, we see that the operator norm (j4.7p in the 
case of the unit disc is precisely \/2. So we see in the case of the unit disc that 
Theorem SH] holds with c = C = ^/2. 

Also note that when I — 1, /i„_i ~ ?i + cn^'^, and then ||wj|| ^ /i . These are 
'whispering gallery modes', which saturate the bound (|4.6p . 

We now sketch some parts of the proof of Theorem 14.11 Let us show how to 
obtain an upper bound for a single function Wj, that is, prove 

This already contains the crucial difficulties in the proof. 

We return to (|2.2[) . and deduce from it, using a vector field V equal to (i„ at the 
boundary, that 

Vjdlvj = 0{nj). 



Ion 

(This is not quite as straightforward as in the Dirichlet case, as one needs to show 
that the left hand side is 0{^^). This requires some integration- by-parts and relies 
on estimate (|4.6I) .) It follows, using {A — fi'j)vj = at dfl, and that A = — d^-|- Ago 
at the boundary, modulo first order operators, that 

Ian 
That is, 

11(1 - h'^Aaa)T^,\\h(dn) - IK^'A^n - l)T^,\\hidn) = 0(1). 
So it remains to show that 
(4.8) Uh'.Asn ~ iff^jWhion) = 0(1) (cf. gID). 

Intuitively, this quantity should be very small, since the operator (/i^Agji ~ 1)+ 
is microsupported in the elliptic region where we expect eigenfunctions should be 
negligible (the wavenumber on the boundary exceeds iij hence waves are evanescent 
in the normal direction). However, there is a difficulty since the microsupport of 
{hjAgn — 1)+ meets the boundary of the hyperbolic region {JTyj < 1}. 



NEUMANN EIGENFUNCTIONS AT THE BOUNDARY 9 

To prove (14. 8|) . we break up Wj spectrally into three pieces. Choose a smooth 
function of a real variable such that (f>{t) = for f < 5/4 and (j){t) = 1 for t > 7/4. 
Then we decompose 

valid for hj sufficiently small, e.g. ft, ■ < 1/2, and therefore 

^■ 

■=1 + 11 + III. 

Roughly speaking, here the first piece / is supported in the frequency range 1 < 
hi < 1 + Cft^/'^, the second piece // is supported in the frequency range 1 + ch^^^ < 
|?7| < 2 and the third piece /// is supported where \r]\ > 3/2. We now estimate 
each piece /, // and /// separately. 

Estimating I. We can estimate this piece purely using L^ spectral theory. To 
do this we observe that on the support of (1 — (j)){h~'^^^ {h'^ Agn — 1)), we have 
{h^Agn — 1)_,_' < 2h}l^ . The required estimate now follows from this and (|4.6p . 

The other estimates are more intricate, and rely on expressing the boundary value 
of Neumann eigenfunctions in terms of themselves and the semiclassical double layer 
potential. More precisely, let Dh denote the integral operator 

Dh{x,y) = -dnyGh{y,x), x^y, 

where Gh{x,y) is the Helmholtz Green function (A — {h^^ + iO)'^)^^{x,y) on K". 
It is well-known that 

Iterating this we find that 

(4.10) w,=D^^Wj, TV = 1,2,... 

According to [H],-D^ is the sum of a semiclassical FIO microsupported in the hyper- 
bolic region {|?7| < 1} in both variables; a pseudodifferential operator of order —A'', 
in the sense that it maps L^(i90) to H^ (dil) with norm 0{h^); and an operator 
microsupported close to {\ri\ = 1} in both variables. 

Estimating III. We use (|4.10p with N = 2, and write /// as 

(ft^Aao - l)+/'0(/if Aao - l)wj = (/i^Aao - l)+/V(^'Aaa - l)Dlwj. 

When we compose (ft^Agji — 1)+ 4>{h^AoQ — 1) with Z?^, we obtain a pseudodiffer- 
ential operator of order —1, that is, mapping L^{dfi) to H^{dfl) with norm 0{h), 
since 4'{h?AgQ, — 1) is microsupported away from {|?7| < 1}, and the estimate again 
follows from this and (j4.6p . 

Estimating II. This is the most delicate estimate, since the spectral cutoff 
(/)(/i|Aao — 1/^j ) is only locahzing frequencies at a distance h'^^^ away from the 
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hyperbolic set {\ri\ < 1} where Dh is an order zero operator. To deal with this we 
use the following lemma: 

Lemma 4.3. The operator 

(4.11) 0(^^!^^)(l-0)(/.^Aa.-l) 
can be represented as an oscillatory integral 

(4.12) (27r/i)-("-i) I e'^^y^y''^'^'^b{y^-^,^,h)d^, 
where 

n-l , 

(4.13) Hy.v'.O = E ^Jk{^^){v-y')k^, 

and h satisfies estimates 

(4.14) \d^dlb{y,th)\<C^0h-^\P\'\ 

The proof will be given in a future article. 

Remark AA. It does not seem possible to write the operator (|4.11l) in the usual 
pseudodifferential form, with phase function {y — y') ■ ^, because then the principal 
symbol would be 

This function loses a factor h~'^'^ when differentiating in either y or ^ and such 
a symbol class does not lead to a sensible calculus (in the sense of having a com- 
position formula, etc). By contrast, with a judicious choice of the Ojk function in 
(|4.13p , one can arrange that the principal symbol of (|4.12p is 

which satisfies the better estimates in (|4.14p , in that there is no loss of powers of h 
when differentiating in y. 

We can write down an oscillatory integral representation for the operator Dh 
as an intersecting Lagrangian distribution. Using this and the oscillatory integral 
representation for operator (|4.1ip given by the lemma, and using (I4.10p with iV = 1 , 
we can write the operator in // as an oscillatory integral involving one factor of 
Dh- In this integral, the phase is non-stationary on the support of the symbol. 
Using integration by parts in a standard way, we can show that this operator has 
an operator norm bound of 0{h}'^)^ which combined with (|4.6p gives the result. 

Remark 4.5. By using (|4.10p with large values of TV, we can show that the L- norms 
oil I &nd III areO(/i°°). 
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Figure 2. Graphs of minimum tp [v] achievable at each E for 
the Neumann case, for il the unit disc, and fow energy, a) F = Id, 
h) F = Ffj^ given by ()5.2p . As in Fig. [l] w is restricted to he in a 
sufficiently large numerical subspace. The curves (shown lighter) 
lying above the lowest show higher generalized eigenvalues of a 
matrix pair used to compute tp (the lowest gives tp itself). 
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Figure 3. Graphs of minimum ip^[v] achievable at each E for 
the Neumann case, for fl the domain shown in Fig. 21 at higher 
energy. Cases a) and b) are as in Fig. [H 



5. Neumann inclusion bound and numerical demonstration 



Using Theorem 14. II as a crucial tool, we propose the following method of partic- 
ular solutions (MPS) for finding Neumann eigenpairs: at each energy E — ji^ we 
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Figure 4. High- lying approximate Neumann eigenfunction of a 
smooth domain V, (density plot shows |i'(x)p), computed by the 
new MPS proposed in Section [S] v minimizes tp [v] at energy 
E — 2096.240170 which is close to a local tension minimum. The 
tension here was tp < 10~^. The eigenvalue is near the 500th. 
Numerical computation used a basis of 400 fundamental solutions 
lying outside f2 (see MPSpack manual [4]). 



minimize the quantity 

(5.1) 

where (cf. Moral (|I3|) 

(5.2) F^{a) 



iFM = 



\Ff,iAsn)id„v)\\L2(dQ) 



the invertible boundary operator is -F),(Agsi) and 



,,1/3 



c < M M 



4/3 



a > n — II 



4/3 



The effect of F^ is roughly to boost the amplitudes of spatial frequencies on the 
boundary which are close in magnitude to the overall wavenumber /x; however, it 
is regularized to limit this boost to a finite value taking heed of the scaling (|4.6p . 
This leads to the identity, analogous to (|3.ip and <\AA\ . 



(5.3) 



(mintF,,H)" 



F^(Aan)-'w;,(F^(Aan)- 



j 



(m^ - lAf 



where again, the min is to be taken over Helmholtz solutions v. First we give a few 
words about why (j5.3|) holds. As before, {xmwip [v\)~ is the operator norm of the 

composite function g i-> / i— >• v, where / = F^{lS.Qii)~^g and v is the Helmholtz 



NEUMANN EIGENFUNCTIONS AT THE BOUNDARY 13 

solution with dnV — f. In a similar fashion to the derivation of p.3p . we have 
^ ^y^ {f,Wj)vj ^Y^ {Ff,{Aon)'^g,Wj)vj 

^y^ {9,Fi,{Asn)''^Wj)vj 

Then a T*T argument, analogous to that in the proof of Theorem 12.11 gives (|5.3p . 



Finally, since F^(Agn)^^ is essentially (1 — /i^Aao)/ , we can use Theorem 14. II 
(together with (14. 6p ) to prove the following tight Neumann inclusion bound (anal- 
ogous to Theorem 13. ip : 

Theorem 5.1. There exist constants c,C depending only on f2 such that the fol- 
lowing holds. Let v he a nonzero solution of (A — /^^)w — in C°^{Vl). Let tp [v] 
be as in (|5.ip . and let Wmin be the Helmholtz solution minimizing tp [v\. Then 

ciF„ bmin] < d(M^,specjv) < Ctp^[v\. 

We postpone the full proof to a future publication. Comparing to the Dirichlet 

case, we note that there are no factors of y E in the bound (i.e. a = 0); this is to 
be expected dimensionally since the Neumann tension tp [v] already contains an 
extra derivative compared to the Dirichlet tension t[u\. 

We end with some preliminary numerical demonstrations of our Neumann MPS 
in 71 = 2 dimensions. The lowest curves shown in a) and b) of Fig. [2] are the 
Neumann analogs of Fig. [1] for the unit disc, comparing two choices of the operator 
F. It is clear that the naive choice i^ = Id leads to large variations in slopes, whereas 
choosing F = F^ given by (j5.2p causes these slopes to become very similar. As 
discussed, the disc allows Neumann modes whose value L^ norms on the boundary 
vary as widely as is possible. Fig. [3] shows, at higher energies, the same but for the 
smooth planar domain shown in Fig. 14] the difference in slopes of the lowest curve is 
less striking. However, we have also plotted curves showing the higher generalized 
eigenvalues relevant to the numerical implementation of the MPS ^. It is clear 
that our proposed F operator causes these higher curves to acquire not only very 
uniform slopes but much less 'interaction' between the curves, both of which should 
lead to an improved numerical method. 

Finally, in Fig.|4]we plot a Neumann eigenfunction computed with our proposed 
MPS using F as in (|5.2p . The tension ip [v] was found to be less than 10~^. The 
constant C in Theorem 15.11 is unknown, but the local slope of tension graph was 
measured to be about 0.5, corresponding to C w 2. Thus the inclusion bounds 
on the eigenvalue are [2096.240168,2096.240172], i.e. about 9 digits of relative 
accuracy. Computation took a few seconds on a laptop, using a basis set of size 
400, and 450 quadrature points on d^. The F operator was approximated to 
spectral accuracy using trigonometric polynomials on the boundary. We note that 
applying F in higher dimensions n > 3 will prove more of a challenge. 
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